load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


;#HERE YOU WILL PUT THE PATH
location="/Users/mypath"

filsw= systemfunc ("ls "+location+"lowice/siconc_allmon_nov*.nc")
filsc= systemfunc ("ls "+location+"control/siconc_allmon_nov*.nc")

fw    = addfiles (filsw, "r")
fc    = addfiles (filsc, "r")

ListSetType (fw, "join")
ListSetType (fc, "join")
print(fw)
print(fc)

;*******************************************;
;Read in variables
;*******************************************
lat2d=fc[0]->latitude
lon2d=fc[0]->longitude

ens =10
steps=20   ; 20 years of November data
tempcss =  new((/ens*steps,292,362/),float)
tempc =  new((/ens,steps,292,362/),float)

do i=0,ens-1

time=fc[i] ->time
ntime= dimsizes(time)

tempsc =  new((/ntime,292,362/),float)
tempsc = fc[i]->siconc

last=ntime-1
first=ntime-20
tempc(i,:,:,:)=tempsc(first:last,:,:)

do m=0,steps-1
im=(i*steps)+m
tempcss(im,:,:)=tempc(i,m,:,:)
end do

delete(time)
delete(ntime)
delete(last)
delete(first)
delete(tempsc)
end do

tempwss =  new((/ens*steps,292,362/),float)
tempw =  new((/ens,steps,292,362/),float)
do i=0,ens-1

time=fw[i] ->time
ntime= dimsizes(time)
tempsw =  new((/ntime,292,362/),float)
tempsw = fw[i]->siconc

last=ntime-1
first=ntime-20
tempw(i,:,:,:)=tempsw(first:last,:,:)

do m=0,steps-1
im=(i*steps)+m
tempwss(im,:,:)=tempw(i,m,:,:)
end do

delete(time)
delete(ntime)
delete(last)
delete(first)
delete(tempsw)
end do

tempcss!0 = "alltimesteps"
tempwss!0 = "alltimesteps"
;************************************************
; average over the simulations/years
;************************************************

  finfo_tcs = tempcss(0,:,:)
  finfo_tcs = dim_avg(tempcss(y|:,x|:,alltimesteps|:))     ; reorder but do it only once [temporary]
  finfo_tws = tempwss(0,:,:)
  finfo_tws = dim_avg(tempwss(y|:,x|:,alltimesteps|:))

;************************
; calculate the anomaly:
;************************

  finfo_ano=finfo_tws
  finfo_ano=finfo_tws -finfo_tcs

  delete(finfo_tws)
  delete(finfo_tcs)
  finfo_ano!0 = "lat2d"
  finfo_ano@lat2d = fc[0]->latitude

  finfo_ano!1 = "lon2d"
  finfo_ano@lon2d =fc[0]->longitude

;***********************************
; perform t-test
;***********************************

  xtmp = tempcss(y|:,x|:,alltimesteps|:)     ; reorder but do it only once [temporary]
  ytmp = tempwss(y|:,x|:,alltimesteps|:)

  delete(tempwss)
  delete(tempcss)

  xAve = dim_avg (xtmp)              ; calculate means at each grid point
  yAve = dim_avg (ytmp)
  xVar = dim_variance (xtmp)         ; calculate variances
  yVar = dim_variance (ytmp)
  sigr = 0.05                        ; critical sig lvl for r
  xEqv = equiv_sample_size (xtmp, sigr,0)
  yEqv = equiv_sample_size (ytmp, sigr,0)
;
   xN=xEqv
   yN=yEqv
   iflag= True                         ; population variance not similar
  prob = ttest(xAve,xVar,xN, yAve,yVar,yN, iflag, False)
  prob!0 = "lat2d"
  prob!1 = "lon2d"
  prob@lat2d =fc[0]->latitude
  prob@lon2d =fc[0]->longitude

  do i=0,291
  do j=0,361
  if( .not.ismissing(finfo_ano(i,j)) )
  if (abs(finfo_ano(i,j)).lt.0.00001)
  finfo_ano(i,j)=finfo_ano@_FillValue
  end if 
  end if
  end do
  end do

;-------------------------
; making the plots
;-------------------------
wks = gsn_open_wks("pdf","NHallanoice_95per_nov_medpolar")

gsn_define_colormap(wks,"BlueDarkOrange18")

  plot = new(1,graphic)                          ; create a plot array

  res          = True
  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res@cnInfoLabelOn = False                     ; turn off cn info label
  res@cnLinesOn           = False
  res@cnFillOn            = True          ; turn on color
  res@lbLabelBarOn        = True           ; turn off individual cb's

  res@mpPerimOn = True
  res@cnLevelSelectionMode =  "ExplicitLevels"

 res@cnLineLabelsOn=False
 res@mpCenterLonF     = 300.

 res@mpFillOn = True
 res@mpLandFillColor = "grey"
 res@mpFillDrawOrder = "PostDraw"

 res@gsnPolar   = "NH"
 res@mpMinLatF            = 50.

  res1          = True
  res1@gsnDraw  = False                          ; don't draw
  res1@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res1@cnInfoLabelOn = False                     ; turn off cn info label
  res1@cnLinesOn           = False; True
  res1@cnFillOn            = True          ; turn on color
  res1@lbLabelBarOn        = False           ; turn off individual cb's

  res@cnMissingValFillColor   = "white"
  res1@cnMonoFillPattern = False ;True
  res1@cnFillPatterns      = (/-1,8,-1/) ;17, 8, 13
  res1@cnLevelSelectionMode =  "ExplicitLevels"
  res1@cnLevels =(/0., 0.05/)
  res1@cnMonoFillColor     = True         ; Use same fill color

  res1@cnLineLabelsOn=False
  res1@cnMonoLineColor= True
  res1@cnLineColor      = (/"grey36"/) ; (/"white"/)

  res@mpGeophysicalLineThicknessF = 2.
  res@lbAutoManage = False
  res@lbLabelFontHeightF  = 0.017
  res@pmLabelBarHeightF  = 0.08


  res1@gsnCenterString = ""
  res1@gsnLeftString = ""
  res1@gsnRightString = ""

  res@gsnLeftString = ""
  res@gsnCenterString = ""

  res@cnFillColors = (/2,3,4,5,7,9,11,12,13,14,15,16/)
  res@cnLevels =(/-60.,-40.,-20.,-10.,-5.,-2.5,0,2.5, 5./)

  res@gsnRightString = ""

plot(0) = gsn_csm_contour_map(wks,gsn_add_cyclic_point(finfo_ano(:,:)),res)
plot1 = gsn_csm_contour(wks,gsn_add_cyclic_point(prob),res1)
overlay(plot(0),plot1)

name1 ="siconc_ano.nc"
name2 ="siconc_significance_ttest95pct.nc"
outf1 = addfile(name1,"c")
outf1->lat2d  = lat2d
outf1->lon2d  = lon2d
outf1->prect  = finfo_ano(:,:)

outf2 = addfile(name2,"c")
outf2->lat2d  = lat2d
outf2->lon2d  = lon2d
outf2->prob   = prob(:,:)
;________________________
; creating the panel:
;________________________
 resP            = True

resP@gsnPanelLabelBar    = False;True                ; add common colorbar
resP@lbLabelFontHeightF  = 0.007               ; make labels smaller

resP@gsnPanelBottom   = 0.05                   ; add space at bottom
;resP@gsnPanelFigureStrings= (/"a)","b)","c)"/) ; add strings to panel
res@txFontHeightF     = .24


resP@gsnMaximize    = True                ; maximize plot
gsn_panel(wks,plot,(/1,1/),resP)


